ANALISI GEOSTATISTICA DEI DATI AMBIENTALI – Livello base (T28_25_3_1)

Francesco Pirotti

2025-09-16

Riferimenti utili

CRS - coordinate reference system
https://r.geocompx.org/spatial-class#crs-intro
https://r-spatial.org/book/02-Spaces.html#sec-crs

meuse dataset
https://cloud.r-project.org/web/packages/gstat/vignettes/gstat.pdf

metodi di interpolazione geografica
https://r-spatial.org/book/12-Interpolation.html

Obiettivo

Stimare i valori di una variabile nel dominio dello spazio (eventualmente anche del tempo)

Cosa si intende per variabilità spaziale?

Dati di esempio

Useremo il dataset “meuse” e le librerie che vedete sotto

library(tidyverse) 
library(stats)
library(ggplot2)
library(terra)
library(gstat) 
library(kableExtra)
library(sp)
data(meuse)
class(meuse)
## [1] "data.frame"
kable(head(meuse, 5), caption = "Dataset *meuse*") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))
Dataset meuse
x y cadmium copper lead zinc elev dist om ffreq soil lime landuse dist.m
181072 333611 11.7 85 299 1022 7.909 0.0013580 13.6 1 1 1 Ah 50
181025 333558 8.6 81 277 1141 6.983 0.0122243 14.0 1 1 1 Ah 30
181165 333537 6.5 68 199 640 7.800 0.1030290 13.0 1 1 1 Ah 150
181298 333484 2.6 81 116 257 7.655 0.1900940 8.0 1 2 0 Ga 270
181307 333330 2.8 48 117 269 7.480 0.2770900 8.7 1 2 0 Ah 380

Cosa facciamo

Vogliamo un modello per “spiegare” la variabilità e dunque poter stimare (interpolare) al meglio dove non abbiamo informazioni.

Stimare valori di una variabile applicando un modello creato conoscendo valori misurati (e.g. campionamenti). E’ ragionevole adottare un criterio spaziale, ovvero usare un criterio di spazio (il valore del punto sarà probabilmente più simile ad una misura vicina che una misura più lontana, a parità di altre condizioni). Da qui la geostatistica.

fonte: Appunti GIS Prof. Alberto Guarnieri - CIRGEO
fonte: Appunti GIS Prof. Alberto Guarnieri - CIRGEO

Giorno 1

Richiamo concetti statistici fondamentali: varianza, covarianza, correlazione, autocorrelazione

Varianza

Vediamo come valutare alcuni momenti di una distribuzione - media e varianza ci danno il primo e secondo momento. Z è la variabile aleatoria, nel caso del dataset MEUSE possiamo scegliere la concentrazione di Zinco.

\[M(Z)=\mu_Z=\mathbb E[Z]\]

\[M(Z)=\mu_Z=\frac{1}{n}\sum_{j=1}^n{Z_j}\]

\[Var(Z)=\sigma^2_Z=\mathbb{E}\Big[\big(Z-\mathbb{E}[Z]\big)^2\Big].\]

\[Var(Z)=\sigma^2_Z=\frac{\sum_{j=1}^n(Z_j-M_Z)^2}{n}\]

Formula più pratica

\[Var(Z)=\sigma^2_Z=\mathbb{E}[Z^2]-\mathbb{E}[Z]^2\]

\[Var(Z)=\sigma^2_Z=M(Z^2)-M(Z)^2 \]

Varianza esempio

Calcoliamo la varianza con tre formule diverse:

MZ<-mean(meuse$zinc) 
VarZ <- var(meuse$zinc)
## media quadratica dei quadrati
MZ2<-mean(meuse$zinc^2) 
VarZ <- MZ2 - MZ^2
VarZ <- mean((meuse$zinc-MZ)^2)
  
 

plot( meuse$y, meuse$zinc, xlab="Coord Y", ylab="Zinco", pch="+")
abline(a = MZ,0 )
rect(min(meuse$y), MZ-VarZ^0.5, max(meuse$y), 
     MZ+VarZ^0.5, col = "#ff000045", border = "#ffffff00")

# varianza  = (dev std)^2

Varianza esempio (2)

hist(meuse$zinc)
rug(meuse$zinc)
rect(MZ-VarZ^0.5, 10000, 
     MZ+VarZ^0.5, 0, col = "#ff000045", border = "#ffffff00")
abline(v=MZ, col="red", lwd=2)

Covarianza

La covarianza tra due variabili casuali \(Z_1\) e \(Z_2\) misura quanto esse tendono a variare insieme. Se positiva variano entrambe nella stessa direzione, se negativa in direzioni opposte.

\[\mathrm{Cov}(Z1,Z2)=\mathbb{E}\Big[\big(Z1-\mathbb{E}[Z1]\big)(Z2-\mathbb{E}[Z2]\big)\Big]\]

Esempio Covarianza

La covarianza dipende dalle unità di misura delle variabili → non è standardizzata. Per questo si usa spesso la correlazione, che normalizza la covarianza tra -1 e 1.

La correlazione = 1 equivale ad una covarianza che ha valore uguale alla varianza.

CovZ<-cov(meuse$zinc, meuse$zinc)

var(meuse$zinc) == CovZ
## [1] TRUE

Covarianza

Analizziamo un gradiente di temperatura simulato

load("oggetti.rda")
plot(lat_seq, temps, type = "l", col = "blue", lwd = 2,
     xlab = "Latitude (°N)", ylab = "Temperature (°C)",
     main = "Gradiente Temperatura: Italia → Danimarca")
points(lat_seq, temps, pch = 16, col = rgb(0,0,1,0.3))

Covarianza

La covarianza equivale alla varianza se la calcoliamo confrontandola a se stessa, ma se trasliamo di un dato intervallo (“lag”) e la ricalcoliamo, la covarianza non diventa zero, diminuisce gradualmente….

cov(temps,temps)
## [1] 6.150969
cov(temps[2:200],temps[1:199])
## [1] 5.122209
cov(temps[3:200],temps[1:198])
## [1] 4.972588
## Provate con una variabile senza correlazione 
# rand <- runif(200) 
# cov(rand,rand)
# cov(rand[2:200],rand[1:199])
# cov(rand[3:200],rand[1:198])

Questo è un concetto molto importante da comprendere.

Covarianza

Più aumentiamo la distanza del lag, come è ragionevole pensare per dati spazialmente correlati, più diminuisce il valore di covarianza (si allontana dal valore della varianza == covarianza ).

Se si abbassa per poi rialzarsi, può indicare un pattern ricorrente (provate sostituendo i valori di temps con un valore ciclico come temps <- sin((1:200)/10) .

nm <- length(temps)
plot(temps[1:200], temps[1:200], pch="x", cex=0.5, xlab="Temperatura", ylab="Temperatura" )
points(temps[1:199], temps[2:200], pch="x", cex=0.5, col="red" )
points(temps[1:198], temps[3:200], pch="x", cex=0.5, col="green" )

Covarianza

Non ci crediamo? Rifacciamolo per un dato casuale con media e varianza uguale a quella del nostro dato:

t.rand <- rnorm(nm, mean = MZ, sd = VarZ^0.5)

plot(t.rand[1:(nm-1)], t.rand[2:nm], pch="x", cex=0.5, xlab="km", ylab="Temperatura" )
cc=1
covs.rand<-list("0"=cov(t.rand, t.rand), "1"=cov(t.rand[1:(nm-1)], t.rand[2:nm] ) )

for(lag in seq(6, 36, 10)){ 
  covs.rand[[as.character(lag)]] <- cov(t.rand[1:(nm-lag)], t.rand[(lag+1):nm])
  points(t.rand[1:(nm-lag)], t.rand[(lag+1):nm], cex=0.5, pch="x", col=cc )
  cc=cc+1
}

plot(names(covs.rand),  unlist(covs.rand), type="b", 
     ylab="Covarianza", 
     xlab="lag" )

Correlazione

Relazione lineare tra due variabili, cosiddetta di Pearson. Semplicemente normalizza la covarianza per le deviazioni standard delle due variabili moltiplicate tra loro:

\[-1\le\rho_{Z1Z2}=\frac{Cov(Z1Z2)}{\sqrt{Var(Z1)}\cdot\sqrt{Var(Z2)}}=\frac{\sum_{i=1}^n(Z1-{\mu_{Z1}})(Z2-{\mu_{Z2}})}{\sqrt{\sum_{i=1}^n(y_i-{\mu_{Z1}})^2}\sqrt{\sum_{i=1}^n(y_i-{\mu_{Z2}})^2}}\le+1\]

Autocorrelazione

L’autocorrelazione misura quanto una variabile è correlata con sé stessa quando viene considerato uno spostamento nel dominio dei valori.

In una sola dimensione — sia essa tempo o spazio (ad esempio nel nostro transetto di temperatura SUD → NORD) — questo spostamento viene detto lag, traducibile come sfasamento.

L’analisi dell’autocorrelazione si effettua tramite la funzione CCF (Cross-Correlation Function), che nel caso specifico valuta la correlazione della serie con sé stessa a diversi lag.

acf() in R ti mostra l’autocorrelazione ai vari lag (sfasamenti).

#ccf(temps, temps, plot = F)
#ccf(temps, temps)

# Autocorrelazione
acf(temps, main = "Autocorrelazione della temperatura")

Richiamo concetti base su sistemi di riferimento coordinate (CRS), geografiche e proiettate

Cosa sono i dati spaziali?

Rappresentare informazioni associate alla loro posizione nello spazio

Due modelli: il modello vettoriale ed il modello raster.

Ogni valore è georiferito, ovvero associato ad un punto in un sistema di riferimento, quindi delle coordinate.

L’accuratezza può variare a seconda dello strumento usato per la misura.

L’accuratezza dipende dalla scala, che determina la tolleranza agli errori.

Modello vettoriale

I dati vettoriali hanno tre primitive geometriche di base: punti linee e poligoni. La base è sempre il punto con una coordinate, in 2d o 3d.

By M. W. Toews - Own work, CC BY-SA 4.0
By M. W. Toews - Own work, CC BY-SA 4.0

Ci sono decine di formati digitali per i dati vettoriali.

Per un riferimento vedete la pagina QUI.

Ad oggi gli strumenti in ambiente R di riferimento per leggere dati vettoriali è la libreria sf.

Modello raster

Il formato raster è una semplificazione/discretizzazione dello spazio in celle/pixel/ nodi a passo costante, quindi una “griglia” georiferita mediante uno o più punti. Il passo della griglia è spesso identico nelle due direzioni (ma non necessariamente), e ne determina la risoluzione spaziale.

fonte: appunti GIS Prof. Alberto Guarnieri
fonte: appunti GIS Prof. Alberto Guarnieri

Modello raster

Può essere 3D (voxel)

Ci sono decine di formati di file che rappresentano dati raster. Per un riferimento vedete la pagina QUI.

Ad oggi la libreria di riferimento per leggere questi dati è la libreria terra, e stars

Attenzione alla gestione della memoria in quanto i dati raster possono diventare molto pesanti dal punto di vista di spazio di memoria (1000 pixel x 1000 pixel è un raster “piccolo” ma già fornisce 1e6 elementi.

Sistemi di riferimento

Sono delle convenzioni - si stabilisce un sistema con una origine ed un orientamento degli assi tra loro perpendicolari - Coordinate planimetriche o polari - Sistemi locali o globali

Sistemi di riferimento - geoide

Prima di tutto bisogna definire un riferimento fisico - la forma della superficie definita come “interfaccia solida” è “complicata” (è elastica e subisce l’effetto di svariate forze che la deformano)… e neanche stabile nel tempo!

Si prende come convenzione la forma che avrebbe la superficie equipotenziale (livello del mare) se non ci fossero terre emerse - lo si chiama GEOIDE. E’ la superficie di riferimento per comprendere dove la gravità spinge l’acqua. NB è in corso con il recente piano di telerilevamento nazionale il calcolo di un nuovo nuovo geoide ufficiale nazionale (vedi qui)

Sistemi di riferimento - ellissoide

Il geoide è fondamentale per comprendere fenomeni legati alle quote, ma molto complesso da trattare. Si semplifica con l’ellissoide, una rappresentazione trattabile matematicamente.

Ci sono tanti ellissoidi che sono stati stimati per avvicinarsi il più possibile alla realtà del geoide

Coordinate geografiche

Latitudine e longitudine, coordinate polari.

Coodinate geocentriche

Coordinate cartografiche

Coordinate planimetriche, su un piano. Ma la terra non è piatta!

Soluzione: accettiamo deformazioni ma le teniamo sotto controllo (scala, accuratezza, tolleranza)

Coordinate cartografiche

Si proiettano i punti nella “sfera” su un piano.

Fonte: QGIS documentation
Fonte: QGIS documentation

I sistemi di riferimento (CRS)

Vedi documentazione di QGIS QUI

Vedi NOTA PER IL CORRETTO UTILIZZO DEI SISTEMI GEODETICI

Un esercizio con R

Le coordinate di Prato della Valle: Lat: 45.398557, Long. 11.876460 trasformate nei vari sistemi. TEST - prova a consultare il documento IGM e trasformare la coordinata nel sistema “RDN2008 / Zone 12 (N-E)” e vedere come differisce dal sistema “RDN2008 / Italy zone (N-E)” … sono uguali?

library(sf)
lat<-45.398557
long<-11.876460
# punti in classe SF
pt <- st_sfc(st_point(c(long, lat)), crs = 4326)
# UTM fuso 32
pt_utm <- st_transform(pt, 32632)
# 
pt_ecef <- st_transform(pt, 4978)
# le coordinate nei vari sistemi
geographic = sf::st_coordinates(pt)
planimetric = st_coordinates(pt_utm)
geocentric = st_coordinates(pt_ecef)

Fondamenti statistici dell’analisi spaziale dati e variabili regionalizzate

Variabili regionalizzate

variabili misurate in punti e correlate spazialmente

Statistica descrittiva di base

Vedi slides iniziali

Metodi deterministici di stima per interpolazione spaziale

Metodi deterministici vs stocastici

I metodi deterministici non considerano la componente aleatoria (incerta, casuale) dei dati, ma si basano su assunzioni deterministiche per stimare i valori in punti non campionati - ovvero non considerano la struttura spaziale, nessuna autocorrelazione e non forniscono incertezza.

Stimano il valore Z* in un punto dove non abbiamo una misura usando le misure note vicine, ma come?

\[Z* = F(x,y) = \sum_{i=2}^n w_i Z_i\] dove il peso del punto \(1\) dipende dalla distanza \(h\) di un \(i-esimo\) numero di punti \(n\)

Vicino più prossimo (Nearest Neighbour - NN)

Semplicemente quello più vicino
\[ Z*= Z_i \] Dove \(Z_i\) è il vicino più prossimo e \(Z_*\) il punto da interpolare

Vicino Naturale (Natural Neighbour - NatN)

poligoni di Thiessen danno l’area di influenza.

\[ Z*= Z_i \] dove il peso è indicato da Sibson

\[ Z_i(\mathbf{x})=\frac{A(\mathbf{x}_i)}{A(\mathbf{x})} \]

dove \(A(x)\) è il volume della nuova cella centrata in x e A(xi) è il volume dell’intersezione tra la nuova cella centrata in x e la vecchia cella centrata in xi.

By Markluffel - Own work, CC BY-SA 3.0
By Markluffel - Own work, CC BY-SA 3.0

Inverso pesato sulla distanza (IWD)

\[Z* = F(x,y) = \sum_{i=2}^n w_i Z_i\] dove il peso dipende dalla distanza \(h\)

\[w_i = \frac {h_i^{-p}}{\displaystyle \sum_{j=i}^n h_j^{-p}}\]

e da una potenza \(p\) che va inserita dall’utente e può essere oggetto di ottimizzazione.

By Hassan Hadji CC-BY-SA Wiki
By Hassan Hadji CC-BY-SA Wiki

Esempio IDW

library(gstat)
data(meuse.grid)
data(meuse)
## molto importanti i passaggi sotto
coordinates(meuse) = ~x+y
coordinates(meuse.grid) = ~x+y
gridded(meuse.grid) = TRUE
### 
zinc.idw = idw(zinc~1, meuse, meuse.grid)
## [inverse distance weighted interpolation]
class(zinc.idw)
## [1] "SpatialPixelsDataFrame"
## attr(,"package")
## [1] "sp"
spplot(zinc.idw["var1.pred"], main = "Interpolazione IDW concentrazione Zinco")

Esportiamo su raster GeoTIFF per GIS

La libreria terra riconosce l’oggetto “sp” - posso aprirlo con qualsiasi sofware GIS! Nota: non vengono specificati sistemi di riferimento.

zinc.idw.terra <- terra::rast(zinc.idw)
plot(zinc.idw.terra, main="Mappa Concentrazione Zinco")

# terra::writeRaster(zinc.idw.terra, "zincIDW.tif")

Esportiamo su raster GeoTIFF per GIS

Conoscendo l’area, vicino a Maastricht - Paesi Bassi - applichiamo un SR probabile e verifichiamo aprendo in GIS

terra::crs(zinc.idw.terra) <- "EPSG:28992"
# terra::writeRaster(zinc.idw.terra, "zincIDW.tif", overwrite=T)

Esercizio finale giorno 1

Eseguite la procedura sopra per interpolare la variabile “lead” (piombo). Salvate il GeoTIFF con il sistema di riferimento corretto nella cartella condivisa QUI dentro la cartella “giorno1” nominando il file cognome_nome.tif